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We explore a class of hybrid (piecewise deterministic) systems characterized by a large number of 
individuals inhabiting an environment whose state is described by a set of continuous variables. We 
use analytical and numerical methods from nonequilibrium statistical mechanics to study the influ- 
ence that intrinsic noise has on the qualitative behavior of the system. We discuss the application 
of these concepts to the case of semi-arid ecosystems. Using a system-size expansion we calculate 
the power spectrum of the fluctuations in the system. This predicts the existence of noise-induced 
oscillations. 
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I. INTRODUCTION 

When modelling complex systems, we often have to 
decide between two contrasting approaches: whether 
we focus on the behavior of the individuals that con- 
stitute the system, or on the aggregate behavior in 
the large. Furthermore, there is usually a trade-off 
between realism and tractability: while individual 
based models (IBMs) are versatile enough to incor- 
porate many details of the system of interest, they 
can soon become computationally unmanageable if 
too many details are included. Moreover, with too 
much detail we may end up sacrificing insight into the 
way that the model works. Similarly, while aggregate 
models, usually based on differential equations for a 
continuous density of individuals, can be far more 
tractable and provide powerful insights through the 
use of analytical tools, they are susceptible to miss- 
ing relevant details that can affect even the qualita- 
tive predictions of the model (see also for a 
recent general review). It seems, then, that an inter- 
mediate approach is desirable in order to find a good 
compromise between realism, efficiency, and insight. 

Stochastic hybrid models [Tsl - fisj are a general class 
of models that can be very useful in exploiting the 
benefits of both of these approaches. Of particular in- 
terest to us here is the subclass of so-called piecewise 
deterministic Markov processes (PDMP) [13, Hsj . 
which has recent ly gain ed increased attention in the 
natural sciences |19^f23j . A PDMP models a system 
characterized by both discrete and continuous vari- 
ables, where the former follow a stochastic process 
while the latter are governed by a deterministic (dif- 
ferential) equation. These two dynamics are coupled 
as the dynamical law for each depends on the current 
state of both. Motivated by applications to ecology, 
here we assume that the discrete variables refer to 
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the number of individuals (e.g. plants, animals, etc.) 
that are immersed in an environment whose state is 
in turn characterized by the continuous variables (e.g. 
amount of water, temperature, etc.). 

Recently, Faggionato et al. [l^ [l^l have initiated a 
study of PDMP by using tools from nonequilibrium 
statistical physics. In particular, these authors have 
obtained some formal results concerning the deter- 
ministic asymptotic behavior of the system, as well 
as the fluctuations about it, when increasing the rate 
of the stochastic transitions of the discrete variables. 
Here, we are also interested in the nonequilibrium 
statistical mechanics of PDMP; however, we will be 
concerned with the qualitative differences in behav- 
ior from the deterministic asymptotic limit, induced 
by the fluctuations, especially the existence of noise- 
induced quasi-oscillations in a regime where the de- 
terministic approximation predicts only a stable fixed 
point. Furthermore, in contrast to the work of Fag- 
gionato et al., we do not restrict the discrete variables 
to a finite set, but allow them to take on any non- 
negative integer value. Indeed, the mesoscopic limit 
we consider assumes a relatively large number of in- 
dividuals, in the spirit of the system-size expansion 
developed by van Kampen . 

The techniques we use are general and can be ap- 
plied to any system that could be modelled in a piece- 
wise deterministic fashion. An example could be an 
environment characterized by a set of continuous vari- 
ables which would evolve deterministically were it not 
for the influence of a finite number of individuals that 
inhabit it. Here we will consider an application where 
the environmental variables are water in an ecosystem 
which contains a population of plants. Traditionally, 
such systems have been modelled by a fully determin- 
istic dynamics of continuous densities. However, fol- 
lowing the discussion above, it would seem more nat- 
ural to model plants as discrete entities, rather than 
to define a density of plants. It is also well known that 
spatial interactions play a relevant role in these sys- 
tems, as the emergence of spatial patterns are com- 
mon [25l - l30j . In this present paper, however, we will 
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only consider non-spatial models, implicitly implying 
that their area is sufficiently large that we can ne- 
glect any spatial structure altogether. A large region 
can of course contain a relatively large number of in- 
dividuals, though still finite. As we will see below, 
demographic noise can still be relevant in such a sys- 
tem, and induce behavior qualitatively different from 
that described by deterministic approximations. For 
smaller numbers of plants these effects are expected 
to be even stronger, emphasizing the importance of 
stochastic modelling in this field. 

The paper is organized as follows. In Section |lT] 
we define the general framework and introduce the 
model for semi-arid ecosystems which will serve as 
an application of the ideas and techniques which we 
develop. In Section IIIII we show how to derive the 
master equation that governs the evolution of the sys- 
tem and in Section ITVl we discuss the idea behind the 
system-size expansion which we use in the analysis. 
In Section |V] we obtain the leading term in this ap- 
proximation which leads to the non-spatial version of 
deterministic e quat ions extensively studied in the lit- 
erature [2l,[13:113- The next-to- leading term in the 
expansion, which characterizes the fluctuations in the 
system, is described in Section IVTl Finally, in Section 
IVIII we present our conclusions. Most of the technical 
details are collected in the Appendixes. 

II. MODEL DEFINITIONS 
A. General setup 

The simplest class of hybrid systems have states 
that are described by a pair of variables (n, x), where 
n is discrete and x is continuous. The former would 
typically be the number of individuals (e.g. plants) in 
the system, while the latter would refer to the state of 
the environment (e.g. amount of water) which these 
individuals inhabit. The dynamics of the continuous 
variables x is deterministic if conditioned on the dis- 
crete variables, n. On the other hand, the discrete 
variables follow a stochastic process whose transition 
probabilities depend on the continuous variables x. In 
other words, the continuous variables are governed by 
a deterministic differential equation of the form 

f.F(„,x). (1) 

By contrast, the discrete variables follow some 
stochastic transition rules that describe the processes 
which individuals, denoted by P, undergo. For in- 
stance, in the specific case of a birth-death process, 
these have the form 

P 2P, (2a) 

P 0. (2b) 



Here Ff, and F^ are the birth and death rates, respec- 
tively, and can depend on the environmental vari- 
ables X. In principle, the transition rules ([2]) can de- 
scribe any other kind of processes such as migration 
or growth. However, in this paper we will focus exclu- 
sively on simple birth-death processes. This consti- 
tutes an instance of a piecewise deterministic Markov 
process (PDMP) [T7l-l20j: the reason for the adoption 
of this name will become clear later. 

If the system is in state (n, x) the transition rates 
are then given by 

Tb{n + l\n; x) = nTbix), (3a) 
Td{n - l\n; x) = nVdix). (3b) 

All of this can be generalized to a system with 
D discrete variables and C continuous variables by 
introducing the state variables n = (711, 71,2, . . . , n^i) 
and X = (xi, X2, ■ ■ ■ , xc)- Having specified the state 
variables and the transition rates, we can now go on 
to write down the master equation that governs the 
evolution of the probability function of the system. 
However, before doing so, it is useful to give a specific 
example of a PDMP in order to make these concepts 
more concrete. 



B. Example: A semi-arid ecosystem 

As an example, we will consider a non-spatial piece- 
wise deterministic model of semi-arid ecosystems, de- 
fined in terms of three variables corresponding to, re- 
spectively, the densities of surface water, cr, soil wa- 
ter, w, and number of plants, n. A spatial version 
of this model has been well-studied in the ecologi- 
cal literature [2^ [13, HH , but in the case where 
the number of plants is effectively infinite. In the 
model, rainwater falls onto the surface of the land, 
and then infiltrates into the soil where it is taken 
up by plants. Although rainfall in a semi-arid envi- 
ronment can vary drastically and unpredictably over 
time scales short in comparison with the birth-death 
dynamics of plants, here we focus on the average 
amount of rainfall over a relatively long period of 
time. This is a simplifying assumption that allows 
us to separate the influence of intrinsic demographic 
fluctuations from extrinsic environmental noise (see 
e.g. js^ - lsBI l for investigations on the latter). Indeed, 
this simplification has also been used by ecologists 
|H, [13, [2^ HH in order to avoid further complications 
related to the infiltration of water under high water 
densities (see e.g. [l^). In comparison to these in- 
vestigations, the only sacrifice of realism in the model 
we study here relates to the neglect of spatial inter- 
actions. In contrast, a realistic feature missing in the 
studies mentioned above, and that we explore here, 
relates to the discrete nature of plants and their in- 
trinsic stochastic behavior. 
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Due to its relative scarcity, water is the main re- 
source that drives the dynamics of semi-arid ecosys- 
tems. It seems natural to model the water den- 
sities X = (cr, w) as continuous variables governed 
by an equation which is deterministic when condi- 
tioned on the number of plants. More specifically, 
we will assume that the water dynamics is governed 
by an equation of the form (Q, but with two con- 
tinuous variables. Following [2^, [l^, [H, [3l| we write 
F={F,,F^) where 

F„{n,:>c) ^ R - a{p) a , (4a) 
F<^(n,x) = a(p)cr - /3(w)p - rw, (4b) 

and where p = p,n, p being a parameter that char- 
acterizes the influence that a single individual P has 
on the dynamics of the environment. In Eq. (0]), i? is 
the average rate of rainfall that increases the amount 
of surface water cr, and r is a constant rate that char- 
acterizes the loss of soil water which can be due, for 
instance, to evaporation (see [26L [37( for an investi- 
gation of a fully deterministic model which considers 
also loss of surface water). Finally, 

a{p) = a — — , /3(a;)=6 — — -, (5) 

p + k u! + k 

are saturable rates that describe the infiltration of 
surface water into the soil and the uptake of soil wa- 
ter by plants, respectively. Here a, b, k and Wq are 
constants. It is worth mentioning that the infiltra- 
tion rate is taken to depend on the number of plants 
in the system, this is in line with [2^, [2^, [3l| and 
reflects the fact that vegetation typically increases the 
propensity of surface water being absorbed into the 
soil. 

We will model plants as discrete entities that fol- 
low a stochastic birth-death process. Later on, we 
shall investigate the impact that demographic noise, 
due to the discrete nature of plants, has on the be- 
havior of the system. We will assume that a single 
plant has an average mass m so that the density of 
biomass per unit area is given hy p ~ mn/A, where 
the integer n denotes the number of plants in the 
system, and where A is the system's total area. It 
is here important to recall that we are focusing on a 
well-mixed system and thus we do not consider any 
spatial structure. The impact that plants have on the 
dynamics of water is via the density p (see Eqs. (Pa]) 
and (|4bp ) . The effect of the creation or removal of a 
single plant on the water dynamics is characterized 
by the parameter p = m/A, the minimal amount by 
which the density p can change due to a discrete event 
of the plant dynamics. The birth and death rates for 
the stochastic plant dynamics arc taken to be 

rfc(c^) = c/3(c^), (6a) 
Td-rf, (6b) 

respectively, where d and c arc constants and /3(w) is 
defined in (O. It should be noticed that the death 



rate is taken to be constant, whereas the birth rate 
is assumed to depend on the current amount of soil 
water. This dependence, and that of F in Eq. (|4|) on 
the density of plants, p = pn, is what couples the 
deterministic and stochastic dynamics of water and 
plants, respectively. 

III. MASTER EQUATION 

We now return to the general development of the 
formalism, and assume a single discrete variable and 
a single continuous variable, to avoid cluttering the 
equations with indices. 

Suppose then that, at time t, the system is in state 
{n',x'), and denote by T{n\n';x') the total rate (i.e., 
including births, deaths, and any other processes in 
the model) for the system to make a transition from 
n' to n. The probability that, in a small time interval 
At, there are no transitions is given by 

pUx') = l-AtJ2mn';x'). (7) 

In this case, the system will evolve detcrministically 
according to Eq. ([IJ, so that after a time At it will 
be in the state (n',a;o) where 

xo^x +F{n,x)At, (8) 

up to first order in At. 

Suppose now that there is a single transition from 
n' to n which takes place precisely at At' < At. In 
this case the system will evolve in a more intricate 
fashion, to the state {n',xi) with 

xi = x' + F{n', x') At' + F(?i, x') {At - At') , (9) 

up to first order in At. Notice that the probability 
to have more than one transition in the interval At 
is ©(At^). 

Therefore, the probability of a transition during a 
time At is composed of two kinds of terms, corre- 
sponding to the two possibilities, ([5]) and above. 
Each of these contain a Dirac delta contribution for 
the pieccwisc deterministic evolution of x. Thus the 
probability of the system being in state (n, x) at time 
t + At, given it was in state (n', x') at time t is given 
by 

VAtin, x\n', x') = PAt{x') 5 {x - xo) 5„',„ 
+AtT{n\n'; x') 5 {x - xi) (1 - 5„' ,n) , (10) 

with xq and xi as defined in Eqs. ^ and (jH) respec- 
tively. Using Eq. ([7]) this gives, to first order in At, 

VAtin, x\n', x') ^ 6{x - xq) (5„',„ 
-At T{i\n'; x')d {x - x') S„,^„ 

+At T{n\n'- x') S {x - x') (1 - (5„/,„) , (11) 
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where in the last two terms we have replaced Xq and 
xi respectively by x' , since these terms are aheady of 
order At. We can now use the Chapman-Kohnogorov 
equation together with Eq. (ITT]) to obtain a master 
equation for the evolution of the probabihty at finite 
time t. The last two terms on the right-hand side 
of Eq. (jlip give the standard terms in the master 
equation, but the first term gives a contribution 

J dx'Y,H^-xo)Sn'^nV{n',x',t). (12) 

n' 

It is worth pointing out that xq depends on the in- 
tegration variable x', see Eq. Introducing a test 
function, and integrating by parts, we find that (e.g. 
see [ssj ) this term equals 

Vin, X, t)-At-^ [F{n, x)V{n, x, t)] , (13) 

to first order in At. This then yields the master equa- 
tion for the evolution of the probability distribution 
P{n, X, t) for the system to be in a state (n, x) (l8l - l20l | 

d d 

Ofin, x,t)=- — [F{n, x)V{n, x, t)] 

+ Y,[Tin\n';x)V(n',x,t) (14) 

— T(n'|n; x)V{n, x, t)] . 

Although we have derived the master equation for a 
PDMP with one discrete and one continuous variable, 
the generalization to many variables can immediately 
be seen to be 

d ^ d 

—Tin, X, t) = - ^ — [F,(n, x)7'(n, x, t)] 

+ J2 [T(n|n';x)P(n',x,t) (^^^ 

-T(n'|n;x)P(n,x,t)]. 

The terms in the second sum are of the standard form 
one would expect in a master equation for birth-death 
processes. The first term describes the (piecewise) 
deterministic evolution of x, and is of the form of a 
drift term in a standard Fokker-Planck equation for 
continuous processes. The absence of a term contain- 
ing second derivatives with respect to components of 
X reflects the (piecewise) deterministic nature of the 
evolution of x. The master equation fully specifies the 
dynamics of the stochastic system, but it cannot be 
solved analytically and it is difficult to solve numeri- 
cally. Instead, either single trajectories for the under- 
lying process can be simulated by using an algorithm 
similar to that originally devised by Gillespie |39l . l40l | , 
see Appendix|B]for details, or approximation schemes 
can be applied to the master equation. In the next 
section we discuss an example of such a scheme. 



IV. APPROXIMATE DYNAMICS 

The method we will use to analyze the master equa- 
tion (fTS)) separates out the average behavior from the 
fiuctuations around it. We will use the example of 
the semi-arid ecosystem in Section III Bl to illustrate 
the idea. In general, the fluctuations in the discrete 
variables n, are expected to be of order ^/n, and their 
impact on the deterministic dynamics of order y/JIn, 
where fi = m/A, as defined above, is the minimal 
change in mass density per unit area induced by a 
birth or death event. We will thus introduce the 
change of variables 

fin = p + ^ri, (16) 

where p is the deterministic density introduced in 
Section III Bl To the order that we will be working, 
p{n) = p, where the angle brackets stand for an av- 
erage with probability density function V(n,x.,t) at 
time t. The stochastic deviation from the determin- 
istic result, given by p, is described by the term ^/Jlr/, 
and so the deterministic limit corresponds to p ^ 0. 
The physical meaning of this limit is that the area of 
the system is to be chosen sufficiently large (formally 
infinite) so that it contains a large number of plants. 
The discreteness of the birth-death dynamics is then 
no longer relevant. In addition the system is always 
assumed to be well-mixed so that any possible spatial 
structure can be neglected. 

The intrinsic fluctuations in the discrete variables 
will induce fluctuations in the continuous variables x. 
For this reason we will also carry out the replacement 

x = X + VmI, (17) 

where % = (x) and ^ are the average and fluctuation 
terms in the continuous variables, respectively. 

We would like to stress that, up to now, Eqs. 
and pT|) are nothing else than a suitable change of 
variables. We are here assuming that all the stochas- 
tic variation is contained in the variables (77, ^), while 
the corresponding averages (p, x), obtained in the 
limit /i — > as explained above, follow a deterministic 
dynamics. For this reason we introduce the probabil- 
ity distribution 11(77, ^,t) of the stochastic variables 
(77,^) alone, in order to separate these two kinds of 
contributions. Once again we reiterate that the fluc- 
tuations are not externally imposed, but are rather 
intrinsic to the model. Their statistical properties 
emerge from the model and from the approximations 
that we carry out, as we discuss in more detail in Sec- 
tion |Vll Notice that, in terms of the new variables, a 
single stochastic transition, n — > ± 1, corresponds 
to 77 — )- 77 ± When p is small, the effect of a single 
transition can, therefore, be conveniently represented 
by a Taylor expansion in ^/u. The detailed calcu- 
lations follow the lines of [24| . and are summarized 
briefly in Appendix El 
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Figure 1; (Color online) Phase diagram of the determinis- 
tic system, Eqs. (|18|l . In the upper region ('Vegetation') 
the state with vegetation p = p* > is stable under 
small perturbations, i.e. Re(Amax) < 0. In the small 
region to the centre left ('Cycles'), there are no stable 
fixed points, but a numerical integration of Eqs. (|18|l re- 
veals the presence of limit cycles (see discussion in Sec- 
tion lVip . In the bottom region ('Desert'), the desert state 
(p = 0) is stable under small perturbations. The two 
points marked 'Simulations' correspond to i? = 1.03 and 
R = 1.05 (with Wo = 0.1), these are the parameters cho- 
sen for the stochastic simulations discussed in SectionjVll 



V. DETERMINISTIC LIMIT 

The leading contribution found from the above ex- 
pansion shows that the average behaviour, described 
by p and x, is given by 



dp 



doc 

dt 



r(p,x) 



where 



(18) 



(19) 



The first equation in ([T5| can be found from Ap- 
pendix 12] or simply by calculating the average of jin 
from the master equation. The function F is given 
by Eq. Q. 

We now analyze some of the properties of the lim- 
iting deterministic dynamics. The first question we 
can ask is: arc there any fixed points, i.e. do states 
exist such that 



dt 



(20) 



and, if so, are they stable under small perturbations 
7 

To study the (linear) stability of these states, it is 
useful to introduce a small perturbation e = {dp, Sx) 
around the fixed point of interest, say (p, x) = 
{p*:X*)-, which is a solution to (^0]) above. We thus 



write (p, x) ~ {p* + ^Pi X* + ^X) and expand equa- 
tions (fT5|) up to first order in e to obtain 



de 
'dt 



J* 



(21) 



Here J* = J{p* ,X*) is the 3x3 Jacobian matrix of 
the system of equations ((T8)) , given by (|A8p , evaluated 
at the fixed point. The stability of the fixed point is 
then governed by the eigenvalue with the largest real 
part, Aniax, of J^* . If the real part of this eigenvalue 
is negative, a small perturbation will die out after 
a time of order l/|Re(Aniax)|, otherwise it will grow 
exponentially fast until non-linearities set in and limit 
the growth. 

We now apply this analysis to Eqs. (jl8p . Depend- 
ing on the choice of parameters, there can be either 
one out of two stable fixed points [l^l or none: the 
fixed points correspond to either a desert state 



po = 0, (To 



R 



R 



a Wo r 



or a state with non-zero vegetation 



cR 
~d 



r ck 
cb-d' 



R 



a{p*) 



dk 
:b-d' 



(22) 



(23) 



Figure [T] illustrates these three situations in terms of 
the parameters R and Wq. In the region where there 
are no stable fixed points a numerical integration of 
Eqs. ([T5|) shows the existence of limit cycles, see Fig.[T] 
and further discussion in Section IVII These limit cy- 
cles were not reported in [2^ , presumably because the 
authors worked in a steady-state approximation that 
neglected perturbations of the surface water density 
fj. As our analysis reveals, these perturbations can 
render the fixed point p = p* unstable. 

While the range of parameters for which there are 
cycles in the deterministic approximation appears 
rather small, it has been observed elsewhere that de- 
mographic noise can effectively enlarge such regions 
[H, [|, II, 0, S IS- Noise- induced oscillations can be 
found in parameter regimes in which a deterministic 
analysis predicts stable fixed points, i.e. outside the 
region labelled 'Cycles' in Fig. [T] These so-called 
quasi-cycles can be characterized by analytical tech- 
niques described below, their amplitude is particu- 
larly pronounced near the deterministic instability. 
Our analysis focuses on parameter regimes outside, 
but near, the region of instability in the determin- 
istic phase diagram of Fig. [T] This is sufficient to 
illustrate the main point we want to make in this 
paper: that demographic stochasticity can alter the 
qualitative behavior of the model relative to that pre- 
dicted in the deterministic approximation. We also 
carry out simulations inside the region of parameters 
in which the deterministic system has limit cycles. 
Finally we would like to add that when spatial inter- 
actions are taken into account, further phases can be 
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Figure 2: (Color online) Power spectrum of the fluctua- 
tions of the number of plants in the phase with a stable 
deterministic fixed point. The solid (red) curve is the an- 
alytical prediction of Eqs. (|A15|) and (|A16|) : symbols show 
results from simulations of the PDMP, averaged over 100 
realizations, with an initial number of plants n = 10''. 
The initial biomass fin, surface water a and soil water 
(jj axe initialized at the fixed point of the deterministic 
approximation. The parameter values used are given in 
Appendix [Cl with R = 1.05 (corresponding to the upper 
point in Fig. [U. 



found in Fig. [TJ In these phases the model exhibits 
spatial patterns coexisting with either of the homo- 
geneous states (desert or homogeneous vegetation re- 
spectively), see m, [13, H^, HH for further details. 
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Figure 3: (Color online) Noise-induced oscillations in the 
piecewise deterministic model for semi-arid ecosystems. 
The continuous (colored) curves show the results of sim- 
ulations with initial number of plants n = 10*; biomass 
H n, surface water a and soil water lj are initialized at the 
fixed point of the deterministic approximation (dashed 
lines). From top to bottom the continuous lines are: sur- 
face water (blue), soil water (red thick line) and biomass 
density (green). The (black) dashed lines show the evo- 
lution of the corresponding deterministic approximation 
initialized slightly away of the fixed point to show the os- 
cillatory convergence. This feature is what couples with 
the noise to induce quasi-cycles in the full model. The pa- 
rameter values used are those given in Appendix [Cl with 
R = 1.05 (upper point in Fig. (TJ. 



VI. FLUCTUATIONS 

Using an expansion in the model parameter fi (ef- 
fectively a system-size or small-noise expansion) it is 
possible to derive a set of coupled linear Langevin 
equations that approximate the stochastic dynamics 
close to a stable fixed point (p*,X*). The expansion 
method goes back to [2J| and is standard by now, see 
[IHE ItI-HiII , so that we do not give full details here. A 
brief summary can be found in Appendix|3 including 
the derivation of the Langevin equation (jAlOp which 
appears in the sub-leading order of the van Kampen 
expansion. The stability of the deterministic fixed 
point is required in order for the stochastic dynamics 
to remain close to the deterministic attractor, so that 
the linear approximation remains valid. 

As detailed in Appendix |21 carrying out a Fourier 
transform of this linear Langevin dynamics with re- 
spect to time allows one to calculate the power spec- 
tra of fluctuations. In particular, we are interested in 

the power spectrum, S{Vl) ~ (\rj{n)\'^'^, of the fluc- 
tuations of the plant density about the deterministic 
fixed point, p* > 0. The quantity rj{n) is the Fourier 
transform (with respect to time) of the fluctuations 
of the plant density, and as before (• • •) denotes an 



average over realizations of the stochastic dynamics. 
Within the van-Kampen expansion analytical results 
can be derived for the power spectrum, S(Jl), de- 
tails can be found in the Appendix, see in particular 
Eqs. (jA15|) and (jA16|) . It is worth pointing out that 
although we here focus on fluctuations of the plant 
density, similar results can be derived for the fluctu- 
ations of the soil and surface water densities. 

We have also carried out numerical simulations us- 
ing a modification of the Gillespie algorithm [sol . l40j , 
as described in Appendix IB] In Fig. [2] we test the an- 
alytical results we have just discussed against direct 
simulations of the PDMP. The value of the parame- 
ter R was set to 1.05. Good agreement between the 
theoretical predictions and the numerical simulations 
is found. The figure demonstrates the existence of 
coherent quasi-cycles, driven by intrinsic noise, in a 
region of parameter space where the deterministic ap- 
proximation predicts a stable fixed point, i.e. where 
one has Re(Amax) < 0. As seen in Fig. [21 the power 
spectrum, S{n), shows a maximum at a characteris- 
tic non-zero frequency, confirming noise-induced os- 
cillations An inspection of a single trajectory, as 
shown in Fig. |31 shows that these quasi-cycles can be 
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Figure 4: (Color online) Effect of demographic noise on 
the limit cycle of the deterministic model. The contin- 
uous curves are results from simulations of the PDMP, 
with initial number of plants n = 10^; biomass ^.n, sur- 
face water a and soil water ui are initialized on the limit 
cycle. Dashed lines represent the deterministic limit cycle. 
From top to bottom the continuous lines are: surface wa- 
ter (blue), soil water (red thick line) and biomass density 
(green). The rainfall parameter is _R = 1.03 (correspond- 
ing to the lower point indicated in Fig. [T}, remaining 
parameter values can be found in Appendix [Cl 



also detected by eye in the time domain. 

In Fig.|4]we show similar simulations for R = 1.03. 
The deterministic system then no longer has a stable 
fixed point, but instead it has a limit cycle. The effect 
of demographic noise on the limit cycle is to introduce 
a stochastic modulation of the period and amplitude 
of the cycle. These effects can be studied analytically 
by separating directions perpendicular to the limit 
cycle from longitudinal modes in a co-moving Frenet 
frame. This is discussed in more detail for a different 
model system in 0, Q • For the present system, how- 
ever, we will not investigate this further analytically, 
but only discuss a few qualitative features. 

Figure[5]shows a projection of the limit cycle (black 
thick continuous line), onto the plane spanned by the 
surface water and soil water variables. In the same 
figure we also plot a part of the stochastic trajectory 
(red thin continuous line) shown in Fig. HI specif- 
ically Fig. [5] shows the segment of the stochastic 
trajectory between time 4200 and time 5500. This 
segment corresponds to the wide plateau of surface 
water concentration cr w 50 seen in the upper panel 
of Fig. m (upper thin solid line). In this segment 
the plant biomass, a measure for the number of in- 
dividuals in the system, remains relatively small, see 
the solid line in the lower panel of Fig. 01 In this 
regime of small numbers of individuals the effects of 
demographic noise are stronger than in periods in 
which there are more individuals present in the sys- 
tem. This is seen in Fig. [SJ the stochastic trajectory 
deviates substantially from the deterministic trajec- 



45 - 




4.5 5 ,5.5 

\ Soil water, CO 




Soil water, co 

Figure 5: (Color online) Projection on the plane of surface 
and soil water variables of the limit cycle (black thick line) 
in Fig.|4l The continuous 'wiggly' (red) curve, near to the 
limit cycle, corresponds to the part of the stochastic tra- 
jectory in Fig. [4] between times 4200 and 5500 (wide valley 
in biomass and wide plateau in surface water). It appears 
that the region where the noise has stronger effects is in 
the top of the plot, which indeed corresponds to biomass 
close to zero (so number of plants relatively small). For 
reference, the dotted (blue) curve at the top corresponds 
to the deterministic approximation initialized at the limit 
cycle, except for the biomass which is slightly displaced 
from the limit cycle (p « 0.003 rather than 0.05). We can 
see that the perturbation initially grows away from the 
cycle and later returns to it. 



tory in the upper part of the limit cycle when the 
plant biomass is small, but it follows the determin- 
istic cycle more closely when the number of plants 
is large (lower part of the limit cycle). As shown in 
[l, Q the overall effect of demographic noise on a limit 
cycle is determined by two different factors: one is 
the relative magnitude of discretization effects, these 
are small for large populations, but more relevant for 
small populations as just discussed. Secondly, the lo- 
cal stability of the deterministic trajectory plays an 
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important role as well. For highly stable deterrain- 
istic trajectories the amplification effect that demo- 
graphic noise undergoes is relatively small. If the de- 
terministic attractor is only weakly stable, then the 
amplification factor can be significant. In the case of 
a fixed point, stability of the deterministic system is 
characterized by the eigenvalue Amax- For limit cycles 
the situation is more complicated, stability is then 
governed by the relevant Floquet exponents 0, II], 
and crucially local stability can vary along the limit 
cycle. In our system we observe that stability in the 
upper part of the limit cycle is only relatively weak. 
To illustrate this we have initialized the deterministic 
system at a point near the limit cycle trajectory, but 
with a slight displacement in the biomass {p ss 0.003 
instead oi p ~ 0.05). No perturbation is applied to 
the water variables a and cj. As shown in Fig. [S] 
(blue dotted line) we find that the perturbation ini- 
tially appears to grow, but ultimately returns back 
to the limit cycle. This relatively weak stability of 
the upper part of the deterministic limit cycle may 
enhance the effects of demographic stochasticity. 



in which a deterministic description would predict a 
stable fixed point. These quasi-cycles can be char- 
acterized analytically by deriving closed-form expres- 
sions for their power spectra within the small-noise 
approximation, in good agreement with simulations. 
This extends existing work on the qualitative fea- 
tures (here quasi-cycles) generated by amplified de- 
mographic noise, see |ll-la.l7l-ll]|. In numerical simula- 
tions we have also investigated the effects of intrinsic 
noise on deterministic limit cycles, similar to the ob- 
servations made in 0, Q we observe that stochasticity 
can lead to a longitudinal phase-drift along the limit 
cycle, coupled with transverse fiuctuations. 

While our present approach is limited to non- 
spatial hybrid systems, work in progress will extend 
these results to piecewise deterministic models with 
spatial interactions. This is particularly interesting 
for more realistic models of semi-arid ecosystems, in 
which pattern forming mechanisms arc known to be 
important [25l - [30| . 



VII. CONCLUSIONS 

In this paper we have studied the effects of demo- 
graphic noise in the framework of piecewise determin- 
istic Markov processes. Specifically we have investi- 
gated a class of hybrid systems, composed of C con- 
tinuous degrees of freedom, and D discrete ones. Such 
systems are of interest for example in the context of 
ecology, where the continuous degrees of freedom can 
represent variables of the external environment such 
as water or light. The discrete degrees of freedom 
could then represent the individual-based dynamics 
of a population of plants or animals. While the pop- 
ulation of individuals follows a standard Markovian 
birth-death process and while the continuous degrees 
of freedom arc governed by ordinary differential equa- 
tions both aspects are not independent. Non-trivial 
behavior arises through the coupling of both types 
of dynamics, for example the birth and death rates 
will depend on the availability of environmental re- 
sources, while these in turn depend on the state of 
the discrete population. 

We have shown how to apply techniques of 
nonequilibrium physics to analyze the effects of de- 
mographic stochasticity on such systems. In particu- 
lar, we have shown how to use a system-size expan- 
sion to obtain a deterministic approximation in the 
limit of infinite populations, and to derive an effec- 
tive Gaussian description of fluctuations about the 
deterministic behavior in the limit of large, but fi- 
nite populations. We have applied these concepts to 
a stylized non-spatial model of semi-arid ecosystems 
and find that demographic noise can induce persistent 
coherent stochastic oscillations in parameter regimes 
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Appendix A: System-size expansion 

In this Appendix wc briefly describe some of the 
intermediate steps of the system-size expansion for 
PDMP in the particular case of the model system 
which we have focused on in this paper. Introducing 
the step-operators [l^l 

e^fin)^ f{n±l), (Al) 

acting on functions /(ti), the master equation (jlSp 
can be written as 

^P(n,x,t) = -Vx- [F(7i,x)7'(n,x,t)] 

+ {e- -l)[n{n + l\n,yi)V{n,^,t)] 

+ (e+-l)[rrf(ji-l|n,x)7'(n,x,0], 

(A2) 

where x = (cr, lj) for our model system. 

In order to make the change of variables described 
in Eqs. (|16p and ((T7)) . wc introduce n(r/,^,t) = 
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V{n,x,t), and use the relation [2J| 

-p{n, X, t) = -n(77, i,t)- — p a^n(r,, I, t) 
1 



^X-V^n(r;,^,t). 



(A3) 



Next the step-operators can be conveniently repre- 
sented by a Taylor expansion in ,JJl 

and similarly, the transition rates can be expanded 
up to first order in ,JJ1. For instance 

M 
A* 

(A5) 

where w stands for the concentration of soil water in 
the deterministic system, and where Tb = Tb(uj) and 

= dT/duj, evaluated at the deterministic value uj. 

Inserting the various ingredients into Eq. (|A2p and 
picking out the term of order p,~^/'^ one finds 

p5„n + x- v^n = $a„n + F- v^n, (A6) 

with $ = x). This equation is satisfied if p and 
X fulfill the deterministic evolution equations (fTS)) . 

Analyzing the next-to-leading term in the expan- 
sion (i.e. the term of order pP) we obtain a Fokker- 
Planck equation 



5tn=|-v^.[j.cn] + ia2[i3n]| 



(A7) 



Here we have used the notation C, ^ (j?,^) 
{il,^a,^Lj)- The quantity 



{dp^p,x) Vx$(p,x) 



(A8) 



is the (3 X 3) Jacobian matrix of the deterministic 
model, and 



B=[d + TbiLo)] p 



(A9) 



is a diffusion coefficient that describes the strength of 
the noise. 

The Fokker-Planck equation (|A7|) is equivalent to 
a linear Langevin equation 



(AlO) 

with a noise term given hy v = {■&, 0, 0) and 

{d{t)i){t')) = B5{t-t'). (All) 



We notice that the second and third components of 
v vanish, all sources of the demographic stochastic- 
ity are now contained in the first component of 
and given by the noise variable ^{t). This reflects 
the fact that only the variable n is subject to demo- 
graphic discretization, whereas x ~ (cr, co) is contin- 
uous. However, given that the three variables of the 
linearized dynamics are coupled (i.e. the Jacobian 
matrix will not generally be diagonal), this does 
not mean that and ^cu will vanish. 

It is also interesting to note that the Jacobian ma- 
trix and the variance B of the noise variable rj 
can be expressed purely in terms of quantities de- 
rived from the deterministic dynamics. During the 
transient phase of the dynamics and B will be 
time-dependent, but in order to proceed analytically 
we will focus on the effect of small fluctuations about 
a fixed point (p, x) = {p* , X*)- The quantities J7 and 
B are then constants. 

Carrying out a Fourier transform of the above 
Langevin equation we find 



with V = (i9,0,0) and 



(A12) 



(A13) 



where the asterisk indicates that the corresponding 
objects have been evaluated at the deterministic fixed 
point. The algebraic equation (jA12p can be solved for 
Cin) to find 



(A14) 



when the inverse exists. This allows us to compute 
the power spectrum of the fluctuations which, by ap- 
plying Cramer's rule, can be written as 



cm 



B 



\det{in-J*) 



\Mism\\ (A15) 



where s = 1,2,3, with ^ = (^, ^cr,$w), and A4is is 
given by the (1, s) minor of the matrix i Ql — J* , i.e. 
the determinant of the matrix obtained after deleting 
its row 1 and column s. Specifically, we are interested 
in the spectral decomposition of the fluctuations in 
the number of plants (i.e. s = 1), for which 



^331 



(A16) 



Appendix B: Simulation algorithm for piecewise 
deterministic Markov processes 

In order to generate trajectories of the piecewise 
deterministic Markov process defined by Eqs. ([T|) and 
([2]), i.e. to sample solutions of the master equation 
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([T5)) , we use a modification of the celebrated Gillespie 
algorithm [s^ . Here we describe one iteration of 
the simulation in the case of a PDMP with one dis- 
crete and one continuous variable. The generalization 
to many discrete and continuous variables is straight- 
forward. Starting at time t with the system in state 
(n,x{t)), to carry out the next simulation step, we 
first draw a time increment r from the distribution 
defined by (see e.g. [ll,[2l|) 

Prob(r < 6') = 1 - cxp ^- T {n; x{s)) ds^ . 

(Bl) 

Here we have written 

T{n;x)^Tbin+l\n;x)+Tdin-l\n;x), (B2) 

which is the total rate for a reaction to occur if the 
system is in state (n, x). The quantity x{s) is given by 
the solution of Eq. ([T]) with initial condition x{t), s > 
t. Generating time increments from the distribution 
above is not entirely straightforward, wc will describe 
the details below. 

Once the increment r has been generated, i.e. the 
time of the next event has been determined, we need 
to decide whether this event is a birth or death event. 
This is carried out following the standard proce- 
dure of the Gillespie algorithm using the appropri- 
ate probabilities given by the reaction rates at time 
t -\- T i.e., a birth event will occur with probability 
Th(n + l\n,Xr)/T{n,Xr), with Xr — x{t -f r), and 
a death event with the complementary probability 
Td(n + ^tItXt) /T(n,Xr)- The event is then carried 
out in the simulation, increasing or decreasing the 
number of plants n by one. So the density is n[t) = n 
prior to the simulation step, and n{t 4- t) = n ± 1 af- 
terwards. Given that the actual birth or death event 
occurs at t+T we have n{t') = n(t) for all t' £ [t, t+r). 
Likewise, the state of the continuous variable at time 
t -f r is given by the solution x{t + r) of the differen- 
tial equation x = F{n,x), with initial condition x{t); 
here the first argument of F is given by the number 
of individuals n at the beginning of the simulation 
step. 

It now remains to explain in more detail how we 
sample the probability distribution defined by Eq. 
(jBip . i.e. how the time increments r are generated. 
This is done by first drawing a number r uniformly 
at random in (0, 1] and then solving 

l-t+T 

J T(n;x(s))ds = ln(l/r) (B3) 



for T. We note, though, that the dependence of the 
left-hand side on the variable r is intricate, so we need 
to use a numerical approximation scheme to evaluate 
the integral while at the same time determining its 
upper limit r self-consistently. Suppose that At is 
a time step suitable for the numerical integration of 
the differential equations x = F(n,x), for every n. 
We introduce the notation xi = x{t -\- 1 At) and Tg = 
T(n,Xi). The algorithm then proceeds as follows: 

1. In the first instance simply approximate the in- 
tegral by Tqt. Substitute this into Eq. (|B3l) 
and determine the resulting value for r. If r is 
found to be smaller than Ai, use the resulting 
value as an approximation for the time incre- 
ment, and carry out the Gillespie step. 

2. If T is found to be greater than At we can im- 
prove the approximation of the integral used 
in step 1: approximate the integral by T^At + 
(t - At)Ti. Substitute this into Eq. (HSI and 
determine r. If r is found to be smaller than 
2At, use this as an approximation for the time 
increment. 



3. Sequentially continue this scheme to the case 
of an increasing number of M — 2,3,... dis- 
cretization steps, i.e. approximate the integral 
by At Efco ^ Tt + {T- MAt)TM and solve Eq. 
(jB3p for T. If r is found to be smaller than 
{M + l)At^ terminate, otherwise increase M by 
one and iterate. 

While this algorithm may appear quite tedious and 
costly in terms of CPU time, we typically find in our 
simulations that the scheme terminates after a few 
steps when /x is small; see also [2l| for a different 
approach. 



Appendix C: Parameter values 

Except when specified otherwise, the parameter 
values used in the simulation are chosen in analogy 
with [li,!!!!!^ and are a = 0.2,6 = 0.05, c==: 10, d = 
0.25, r = 0.2, = 5 and Wo = 0.1. 
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